Solving the boundary value problem for finite Kirchhoff" rods 
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The KirchhofF model describes the statics and dynamics of thin rods within the approximations 
of the linear elasticity theory. In this paper we develop a method, based on a shooting technique, 
to find equilibrium configurations of finite rods subjected to boundary conditions and given load 
parameters. The method consists in making a series of small changes on a trial solution satisfying 
Ph' the Kirchhoff equations but not necessarily the boundary conditions. By linearizing the differential 

equations around the trial solution we are able to push its end point to the desired position, step by 



^ : step. The method is also useful to obtain configurations of rods with fixed end points but different 

• mechanical parameters, such as tension, components of the moment or inhomogeneities. 



PACS numbers; 02.60.Lj, 46.70.Hg, 87.15.He, 87.15.La 



I. INTRODUCTION 



The study of conformations of slender elastic rods is of substantial utility in several applications, ranging from the 
^ ^ fields of structural mechanics and engineering to biochemistry and biology. Examples are the study of coiling and 
Q . loop formation of sub-oceanic cables [|], ||, Q , filamentary structures of biomolecules [|[ ||, 0, || and bacterial 
' fibers [l^j [Tl|], the phenomenon of helix hand reversal in climbing plants and the shape and dynamics of cracking 
whips |l3| . 

• The KirchhofF model |Q provides a powerful approach to study the statics and dynamics of elastic thin rods ^ . 
Ph[ In this model the rod is described by a set of nine partial differential equations (the Kirchhoff equations) in the time 
and arclength of the rod. They contain the forces and torques plus a triad of vectors describing the deformations of 

' the rod. These equations are the result of Newton's second law for the linear and angular momentum applied to the 
^ thin rod plus a linear constitutive relationship between moments and strains. The Kirchhoff model holds true in the 



approximation of small curvatures of the rod, as compared to the radius of the local cross section |16 . An interesting 



1^ , characteristic of this model, known as Kirchhoff kinetic analogy, is that the equations governing the static problem 



are formally equivalent to the Euler equations describing the motion of spinning tops in a gravity field. The Kirchhoff 
equations for equilibrium configurations can, therefore, be written in hamiltonian form. 
' The Kirchhoff equations have been solved for a number of simple situations. Shi and Hearst Q first obtained 
] analytical solutions of the equilibrium equations and, recently, Nizzete and Goriely [0 completed the study by 
making a classification of all kinds of equilibrium solutions. Goriely and Tabor |l^ developed a method to study 
Q the dynamical stability of these solutions and Fonseca and de Aguiar applied this method to study the near 
' ^ , equilibrium dynamics of non-homogeneous closed rods in viscous media. Recently, Tobias et al. |^ developed the 
^i-j' necessary and sufficient criteria for elastic stability of equilibrium configurations of closed rods. 

, In many cases of interest, including biological molecules, the filaments are subject to boundary conditions. Examples 
Ph' are the problem of multiprotein structures, such as histones and gyrase, about which long pieces of DNA wrap [22[, 
^ ' and multiprotein structures, such as the lac repressor complex [^^. Despite the many achievements described above, 
the study of the boundary value problem (BVP) associated with Kirchhoff filaments is still a big challenge. While 
the integration of differential equations from initial conditions is a relatively simple numerical task, the difficulties of 
finding solutions for given boundary conditions are well known in classical mechanics, electromagnetism and quantum 
mechanics. Typical examples are classical trajectories connecting two given space points in the time t, electric 
potentials that vanish at given surfaces and eigenvalues of the Laplace operator defined inside a finite domain (quantum 
billiards) . 

Because of the analogy with spinning tops, the case of trajectories of Hamiltonian systems is of particular interest 
here. The monodromy method, developed by Baranger and Davis p^ , was designed specifically to find periodic 
solutions of Hamiltonian systems with degrees of freedom. Xavier and de Aguiar extended the method to find 
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non-periodic trajectories with any given combination of 2N position and momenta at initial and/or final times. 

A widely used method to solve BVPs is the so called shooting method (26[ . For a single second-order differential 
equation, the method consists in finding the proper 'velocity' at the initial point so as to reach the desired 'target' at the 
end point, similar to the shooting of a projectile. Examples of applications of this method to the Kirchhoff equations 
arc the search for homoclinic orbits in reversible systems hetcroclinic orbits resembling tendril pervesions p9| ] 
and the study of localized buckling modes of thin elastic filaments ^ . 

In the case of open rods some specific BVPs were solved recently. Karolyi and Domokos jS^, using symbolic 
dynamics, found global invariants for BVPs of elastic linkages, as natural discretization of continuous elastic beams, 
an old problem solved by Euler(see reference Gottlieb and Perkins |Q investigated spatially complex forms 

in a BVP governing the equilibrium of slender cables subjected to thrust, torsion and gravity. Also, the criteria of 
Tobias et al. was applied to linear segments subjected to strong anchoring end conditions, where not only the 
end points but also the tangent vector at the end points are held fixed. The dependence of DNA tertiary structure 
on end conditions was studied in [ p2| , where explicit expressions for equilibrium configurations were obtained for a 
specific case with symmetric end conditions. 

Our aim in this paper is to develop a method to find equilibrium solutions of finite rods subjected to boundary 
conditions at their end points and with given load parameters. We emphasize that this is different from the approach 
in [0 where the authors use shooting methods to calculate localized buckling modes. These modes are treated 
as homoclinic solutions of the Kirchhoff equations, corresponding to infinite rods that become asymptotically straight 
in the infinite. Our objective is to find equilibrium solutions for finite rods subjected to boundary conditions at both 
ends. 

Our method is an adaptation of the monodromy matrix method for non-periodic trajectories [ p5[ to the hamiltonian 
formulation of the static Kirchhoff equations. We work with the Kirchhoff equations directly in Euler angles, instead 
of using the Cartesian position and tangent vectors 0. The difficulty in working with Euler angles is that the 
variables that one wants to hold fixed, the spatial position of the filament end points, are not the variables appearing 
in the differential equations describing the rod. However, the number of differential equations to be solved is much 
smaller in these variables. Using a symmetry of the Hamiltonian, we end up with only two independent equations to 
solve. 

One of the motivations of this work is its possible biological applications as, for example, the study of single DNA 
molecules manipulated by optical traps ]3^ , |35| , |3^ , and the DNA loops between multiprotein structures (such as the 
lac repressor-operator complex) ||3^, . 

This work is organized as follows. In Sec. II we review the Kirchhoff equations and, in Sec. Ill, their hamiltonian 
formulation. In Sec. IV we describe our method for solving the BVP. The monodromy method enters as part of 
the solution, and proves to be a very efficient tool. In Sec. V we give numerical examples, calculating the three 
dimensional configuration of rods with different sets of load parameters and end positions. We also discuss the 
existence of solutions as function of the load parameters. In Sec. VI, motivated by the repeated sequences of base- 
pairs commonly found in DNA molecules, we consider rods with periodically varying Young modulus. We compare 
the configurations of these non-homogeneous rods against their homogeneous counterparts, fixing the same end points 
and mechanical parameters. In Sec. VII we summarize our conclusions. 



II. THE KIRCHHOFF EQUATIONS 



The Kirchhoff model describes the dynamics of inextensible thin elastic filaments within the approximation of linear 
elasticity theory They result from the application of Newton's laws of mechanics to a thin rod, and consist of two 
equations describing the balance of linear and angular momentum plus a constitutive relationship of linear elasticity 
theory, relating moments to strains. The Kirchhoff model assumes that the filament is thin and weakly bent (i.e. its 
cross-section radius is much smaller than its length and its curvature at all points). In this approximation it is possible 
to derive a one-dimensional theory where forces and moments are averaged over the cross-sections perpendicular to 
the central axis of the filament. 

A thin tube can be described by a smooth curve x in the 3D space parametrized by the arclength s, and whose 
position depends on the time: x — x(s,t). A local orthonormal basis, (or director basis) d,; = di(s,t), i = 1,2,3, is 
defined at each point of the curve, with chosen as the tangent vector, = x'. In this paper we shall use primes to 
denote differentiation with respect to s and dots to denote differentiation with respect to time. The two orthonormal 
vectors, di and d2, lie in the plane normal to da, for example along the principal axes of the cross section of the rod. 
These vectors are chosen such that {di, d2, d^} form a right-handed orthonormal basis for all values of s and t. The 
space and time evolution of the director basis along the curve are controlled by twist and spin equations 

d- = kxdj, d, =a;xd, i = l,2,3 (1) 
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which follow from the orthonormality relations • dj = Sij . The components of k and u) in the director basis are 
defined as k = X]i=i ^^'^ — Si=i'-^idi. fci and ^2 are the components of the curvature and fca is the twist 
density of the rod. The solution of the twist and spin equations determines 6.3(3, t), which can be integrated to give 
the space curve x(s,i). 

Let the material points on the rod be labeled by 

X{s,t) =K{s,t)+r{s,t), (2) 

where 

r(s,i) = di(s,t) +X2 d2(s,t) (3) 

gives the position of the point on the cross section S{s), perpendicular to x'(s), with respect to the central axis. The 
total force F = F{s,t) and the total moment M = M(s,t) (with respect to the axis of the rod) on the cross section 
are defined by 

F = / psdS. (4) 



M = [ rxpsdS, (5) 

Js{s) 

where Ps is the contact force per unit area exerted on the cross section S{s). In terms of the director basis we write 

F = ELi /^dz and M = ^^^^ M,d,. 

In order to derive a set of equations describing the rod as a one-dimensional object, the rod is divided into thin disks 
of length ds and cross section S{s). To each of these disks the conservation laws of linear and angular momentum are 
applied [0. The result is 

F' + / fe.t dS^ I poX dS, (6) 

Js{s} JS{s} 



M' +x' xF+ [ r X fe^t d5 = / pqv x X dS. (7) 

Js{s) Jsis) 

where fe^t is an external force that will not be considered in our calculations {(ext = in what follows). 

In this article we are interested only in the equilibrium solutions and, therefore, we shall drop the derivatives with 
respect to time. Assuming that the rod has a uniform circular cross section of area A, Eqs. (^) and (|^) can be 
simplified to yield 

F' = 0, (8) 



M' + ds X F = (9) 

which are a set of six equations for 9 variables: F, M and k (from which we determine d^). In order to close the 
system of equations we need a constitutive relation relating the local forces and moments (stresses) to the elastic 
deformations of the body (strains). In the linear theory of elasticity, and for a homogeneous elastic material, the 
stress is proportional to the deformation. The Young's modulus E and the Shear modulus p characterize the elastic 
properties of the material. Therefore, it is possible to obtain, for small deformations, a constitutive relation for the 
moment. For a isotropic material, in the director basis, this relation is JTot : 

M^EI (fci - fcj') di + EI (fc2 - ^2 ) d2 + 2/1/ (fcs - fcs ) da, (10) 

where / is the principal moment of inertia of the cross section, ki are the components of the strain vector and fc" 
are the components of the twist vector in the unstressed configuration. The case fc" = corresponds to a naturally 
straight and untwisted rod. We shall assume /c" — 0. 
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Eqs. (^), and (|T^) can be further simplified by the introduction of scaled variables: 

s ^ sL, F ^ Ffl, 

(11) 

In the new variables the rod has total length L = 1. The static Kirchhoff equations become 

F' = 0, (12) 

M' + d3xF = 0, (13) 

]V[ = fcidi + fc2d2 + rfc3d3. (14) 

where F — 2^/ E is an elastic parameter that does not affect the equilibrium solutions. 

III. HAMILTONIAN FORMULATION 

In order to construct a hamiltonian formulation of the Kirchhoff equations we first note that Eqs. (|l2|)-(|l^) are 
integrable if E and /i are constant iQ. Eq. ( p^ ) shows that the tension F is constant. Let us choose the direction of 
the force as the z direction: 

F = Fez- (15) 
In analogy to the spinning top, the tension F corresponds to the gravity field — mg. Here, F can be considered as an 



external parameter and not as a first integral. Substituting Eq. ( |15[ ) in Eq. (13) and projecting along ez we get 

M' • ez = M'z = (16) 

which does represent a first integral. By projecting the Eq. (|l^) along ds we obtain another integral, M3, since 

M' • dg = M;^ = 0. (17) 

Finally, it is also possible to show that the elastic energy per unit arclength 

i7=iM-k + F-d3 (18) 

is constant, i.e., H' = 0. Therefore H is the last integral. 

The orthonormal Cartesian basis can be connected to the director basis by Euler angles with 

3 

d,-^5.,e, (19) 



where 



cos 6 cos (j) cos ip — sin sin i/j cos 6 cos (j) sin + sin (f) cos i/j — cos sin t 
S~ I — cos 6* sin (/i cos -0 — cos sin ■0 — cos 6* sin sin ■(/; + cos cos ■0 sin sin 6* | . (20) 
sin 6* cos?/; sin sin?/; cos 9 



The static Kirchhoff equations (|l2|)-(|lj) can then be written in terms of 9, (j) and ip. We get 

9" - (0')^ sin 61 cos 61 + F0'(0' + cos 61) sin 61 = i^sin6' 

?/'"sin6' + 2V;'6l'cos6l-F6''(0' + V''cos6') = . (21) 

0" cos 9 ^ijj' 9' sin6'-0" 



These equations can also be derived directly from Eqs. (|16D-(|18D. In terms of the Euler angles the Hamiltonian 
becomes 



where 



Z ZL 2 sm & 



Pe = 0\ (23) 
P^ = M3=r(0' + V/cos0), (24) 

=Mz^ ^' sin^ 9 + cos 6. (25) 

Eqs. (|2^) correspond to Hamilton's equations = and a' — -^p- for a — O^ip oy (f>. We see immediately that 

and P^ are constants and that 9 is the only independent variable. 
The total elastic energy of the rod can be calculated by the integration of the Eq. (|2^): 

Et= ( H{s)ds, (26) 
Jo 

The energy is a function of P^, P^ and P. It also depends on the initial conditions 6{s — Q)=9q and Pe{s = 0) = Peo- 
The procedure to construct equilibrium solutions for given constants P^ and P4, and initial condition (^c-Peo) is 
as follows: first we solve the equations P'g = — §^ and 0' = to obtain {0{s), Pe[s)). Second, using Eq. (|2^), we 
obtain il]{s). The solutions 6'(s) and -0(5) are sufficient to construct the rod by integrating the tangent vector d-^: 

x(s) = / d^{s')ds' . (27) 

Explicitly, 

x{s)^[ am0{s')cos'ip{s')ds' , (28) 

2/(s) = / sin6'(s')sinV'(s')'^s' : (29) 



(30) 







z{s) = / COS 6'(s')ds' , 
Jo 



where 

-R/j — P<t, cos 9(s' 
sin^ 9(s 



V'(s) =^0+1 . ^ ds . (31) 



Substituting equation (^) in (|2^(^o|) and re-arranging the terms we obtain, in matrix form, 

cosV'o sinV^o 0\ / xo{s) 

-sinV^o cosf/'o ?/o(s) I • (32) 
1 / \ z{s) 

where xo{s) and j/o(s) are the equations ( |28| ) and ( p9| ) for ■f/^o = 0. Therefore, it suffices to find the solution with 
ipo = 0. The solutions for other values of ipo are simple rotations of this basic solution. 
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IV. THE LINEARIZED METHOD 

In this section we present our method for finding the configuration of finite rods subject to boundary conditions in 
the position of its initial and final points. In fact, since the Kirchhoff equations are invariant under space translations, 
we can always choose the initial point be the origin. As we saw in the previous section, equilibrium solutions for the 
static Kirchhoff equations depend only on two initial conditions, namely, Oq and Pgo- The third initial condition, ipQ 
corresponds to a rotation of this solution around the z-axis. The problem is then that of finding a solution that starts 



from the origin and ends at zj and at a distance r/ = yyj + xj from the z-axis. Our method is based on a series 

of small deformations made upon an initial solution of the Kirchhoff equations which, however, does not have the 
desired boundary values for rf and Zf. We call this initial solution the trial solution. The method consists in pushing 
the end point of the trial solution to the desired position, step by step. The basic idea, which is a type of shooting 
procedure, is to find a variation in the initial conditions so as to obtain the desired variation in the end point. 

In order to do so, we shall employ a variation of the Monodromy Method originally devised to calculate 

periodic solutions of chaotic Hamiltonian systems. As discussed in the Introduction, working with the Kirchhoff 
equations in Euler angles poses an extra difficulty on the already hard problem of satisfying boundary conditions: the 
variables to be held fixed, r/ and Zf, are not the ones entering the equations of motion, namely, 9, Pg and "tp. The 
advantage is that we can find the solutions solving only two differential equations. 



The rod can be obtained from the Euler angles the equations (|28|), (g9|) and (|30|). The Euler angles, in their turn, 
obey the equations 

e' = Pe (33) 



^, ^ jp, -p, cos e)p, ^ iP,-P,cosefco.e ^ 

sin 6 sin 6 



I p^,-p^cose{s') ^ , 

sm 0[s') 



e 

and 



If we integrate Eqs.(|33|)-(|35|) using the initial condition provided by the trial solution and further integrate Eqs. ^ ^ 

( ^ ) with the resulting Euler angles, we get, of course, the trial rod. Variations in these initial conditions will produce 
variations in the rod configuration, and, in particular, in its end point. In what follows we shall construct an explicit 
relation between a small variation in the initial variables 6q and Pqq and the rod's end point, represented by r/ and 
Zf. Explicitly, we shall find the matrix B such that 

Once B is obtained (and if it can be inverted) we can work our way from the trial solution, whose end point is at, 
say, rt and Zt, to the desired end point at and 2/, provided we do that in a series of small steps. In each step we 
use the previous solution as the trial input, pushing the rod's end point slowly towards its final destination. 



Using r/ — \ x^ + yj, the components of the matrix B can be written as: 



_drf _ Xf dxf Vf dyf 



B,, = -^ = ^pL + yipL (38) 



dPeo rf dPeo rf dPeo 



(40) 



From Eqs. (P§)-(pd) and (BSh we find 



Sxf — / cos 0{s) cos ip{s)S9{s)ds — / sm6{s)sm^p{s)S^p{s)ds, 
Jo Jo 

5yf — I cos9{s)s\\iilj{s)56{s)ds + I si\iO{s) cos%l^{s)5%l^{s)ds, 
Jo Jo 

5zf = -f sm0{s)Se{s)ds, 
Jo 

and 

6i^{s) = / A{9{s'))S9{s')ds' , 
Jo 

where A{9) is given by 

A(e) = 1^- 2(^^-^0Cosg)cosg 

Finally, to find the relation between the variations 59{s) and 5Pe{s) and their values at the initial point s = 0, 
consider small variations of Eqs.(^3|) and around the trial solution : 



59' = 5Pe, ,^ 
6P'^ = C{9)b9, 



where C{9), given by 



2 {Pf~P4,cos9)iP^-AP^cos9) ^{P^^^P^cosflcofl 

C(9) = -P^ r^- r^- h cos9, 

sm 9 sm 9 

is computed at the trial solution. 

The solution to these linear equations can be written in matrix form as 

69(3) \ _ f Mni-s) M,2is) \ ( S9o 



6Pe{s) J \ A'hiis) M22(s) ) \5Po ) ' 

where M is the tangent matrix, satisfying M(0) = 1. In the special case where the trial solution is periodic, M 
called the monodromy matrix. 
Writing 59{s) explicitly as 

59{s) = Mn{s)S9o + Ah2{s)5Po, 
and using Eqs.(|4l|)-(|4^) we can readily obtain 

cos9{s)cosip{s)Mii{s)ds- [ sm9{s)smi;{s) [ A{9{s'))Mii{s)ds'ds 

dxf 



I cos9{s)cos'ip{s)Mi2{s)ds ~ sin 9 (s) sin ip{s) A{9{s'))Mi2{s')ds'ds 
"Peo Jo Jo Jo 



cos9{s)sin^{s)Mii{s)ds + / sin 6'(s) cos -0(s) / A{9{s'))Mii{s')ds' ds 



dOo Jo ^ 

pL ^ f cos9{s)sin^{s)Mi2{s)ds+ [ sin 6'(s) cos V^s) / A{9{s'))Mi2{s')ds'ds 



dPeo Jo 
dzf 



f'^o Jo 
dzf 

dPeo 



1 

sin9{s)Mii{s)ds 



- I sin9{s)Mx2{s)ds 
Jo 
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and, therefore, the matrix B, 

Since we hnearized the equations of motion, we have to check if the new solution, starting from 6i = 6q + 66q and 
Pgi = Pq + SPgQ generates a rod with the chosen final point, within a given precision. If the precision is not reached, 
we can use the newly computed solution as a new trial one, using again Eq. (|3^), now with {Srf, Szf) corresponding 
to the distance between the fixed end point and the end point of the previously computed rod. The process can be 
repeated until the desired accuracy is obtained. 

Finally we note that the elements Mij{s) can be computed by solving the linear equations ( ^ ) with proper initial 
conditions. Indeed, setting 59q = 1 and 5Peo = 0, Eq. ( |48| ) gives Mii(s) = 50{s) and M2i{s) — SPg{s). If, on the 
other hand, we set 69o = and SPgo = 1 we get Mi2(s) = S9{s) and M22(s) = SPe{s). Therefore, Afii(s) and M2i(s) 
are solutions of the linearized Eqs. ( ^6| ) with the initial conditions SOq = 1 and SPgQ = and Mi2(s) and M22(s) are 
the solutions of the same equations with S9q — and SPqq = 1. 

In many cases we might want to push the rod's end-point to a position ry = (r/,z/) far from that of the initial 
trial solution, = (rt, zj). To do that we can divide the Hne connecting to r/ into small segments and apply the 
linearized method N times, moving a small distance at each step. The number of steps required will depend on the 
particular configuration and possibly on the stability of the rod. In all integrations presented in this paper, we used 
a fourth order Runge-Kuta method with fixed step. In all cases the distance between the end-point of the trial rod 
and the target position was divided into 10 segments and the solution converged to the desired boundary condition 
with a precision of 10^^ in each component r and z. 

Figure |l] shows a example of the method. We have chosen the following load parameters in scaled units (see Eq. 
(|ll|): P^ ^ 1.0, P^ = 1.0 and F = 1.0. The desired end point is r/ = 0.7 and z/ = 0.0. We plot the trial, an 
intermediate and the converged rods together, in order to show the process of convergence from the trial to the 
desired solution. The trial solution was computed integrating the Kirchhoff equations starting from 9o — 0.5 and 
PgQ = 1.0, which corresponds to a rod whose final point is ry ~ 0.8 and Zf ~ 0.455. The intermediate solution was 
computed from 9o ~ 0.216 and Pgo ~ 2.257, which corresponds to a rod whose final point is r/ ~ 0.73 and Zf ~ 0.15. 
Finally, the initial conditions obtained for the converged solution are 9o ~ 0.307 and Pgo ~ 2.510. 



V. NUMERICAL EXAMPLES 



The particular trial solution used in the previous numerical example converged smoothly to the chosen final position. 
In some cases, however, a given trial solution does not converge to its destination no matter how many intermediate 
steps are used to divide the line between rj and r/. As we shall see, this problem has to do with the existence or not 

of solutions for a given r/. In scaled units, it is obvious that there are no solutions if + zj > 1. The restrictions 

are actually much stronger than this simple 'length rule', and depend on the values of P^, P^ and F. It might also 
happen that the solution for a given ry does exist, but that the straight line connecting to r/ passes through 
forbidden regions, hindering the convergence. In this section we investigate the space of possible solutions and give 
several examples of rods computed with our method. 

Each initial condition 9q and Pgo leads to an end point r^. The easiest way to map all possible final points is to scan 
the space of initial conditions. Therefore, for a fixed set of parameters P^, P^ and F we calculate Vf = rf{9o,Pgo) 
and plot the resulting figure in the {rf,Zf) space. Points outside this region are unreachable by the rod. Changing the 
parameters, such as the tension, changes the region of possible solutions, including end-points that were not previously 
present and excluding others. 

The results in this section are presented as follows: for each fixed set of P^, and F we show the (r/, z/)-space 

of possible solutions. On the same plot we draw curves of constant D = ^r'j + z'j and, for each D we compute the 

three dimensional configuration of a few rods. 

In order to compare the rods, we always adjust the value of ■00 such that the rod ends va. yf = plane, ipo is 
determined by: 

yO 

tanV-o^-n- (51) 
^/ 

where x'j: and y^ are the final values of x and y for the rod calculated with ■00 — 0. 

In all cases tested our method converged with at least six significant figures to the previously defined final values 
of rf and Zf. 
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A. = 0, = 1 and F = 0.1 

Figure ^ shows the {rf,Zf) map for this case. The map was generated by varying the initial condition in the 
intervals < 6*0 < tt (30 points) and —5.0 < Pgo < 5.0 (200 points). For larger values of Pgo the total elastic energy 
of the rod increases and the final points tend to concentrate in the region near the origin (data not shown). The full 
thick line represents the curve D = 1.0, which is the natural limit for the solutions. But there is a large region inside 
the D = 1 line where no solutions exist. We shall compare it with that of other load parameters later on. We also 
show the lines of constant D for D = 0.9 (full line), D = 0.7 (dashed line) and D = 0.4 (dotted-dashed line). It is also 
interesting to note a forbidden region centered around z^^ ~ and rf 0.25. From the sequence of points crossing 
each other, it is evident that there are two sets of initial conditions {do, Pea) that generate rods with the same end 
point. In general they correspond to rods that are above or below the z axis, and we shall call them the 'up' solution 
and the 'down' solution respectively. 

Figure ^ shows examples of the rods whose end points are marked with circles in Fig. |^. We show the up and 
down solutions for each final point. Figs. |^ and|^b show three rods each for D = 0.4 and D — 0.7. Fig. ||c; shows 5 
different rods for D — 0.9. 

Notice that when or P^ axe zero, the hamiltonian (|2^), becomes symmetric under the transformation F — > —F 
and 9 ^ TT — 9. In these cases the 'up' solution for a given F and is identical to the 'down' solution for —F and 
TT — 9q. When both P^ and P,p are non-zero the symmetry disappears. 

B. = 0, P,^ = 5 and F = 1 

The (rf,Zf) map for these parameters is shown in Fig. ^ It has a curious pattern of thin bulbs centered around 
the z/ = axis, that degenerate for small r/. The only rods possible in this case are those that return close to the 
z = plane. 

Keeping P^ — Q and F = 1 but increasing results in even thiner bulbs. Figures || (a) and (b) show the up and 
down solutions, respectively, for r/ = 0.7, Zf — Q and = 1, 5 and 10. Panels (c) and (d) show the up and down 
solutions for the same parameters, except for rf — 0.9. Notice that large values of P^ correspond to horizontal helical 
rods. 



C. Pv, = P4, and P = 1 

Figure || shows the {rf,Zf) map for P^, = P^ = 5. It resembles the map on Fig. ^, only distorted towards Zf — I. For 
these values of the parameters the rod admits 'vertical' configurations, as opposed to the 'horizontal' configurations 
displayed in the previous case. Figures ^ show examples of rods with P^ — P^ — h and P^ = P^ = 10 for Vf = 0.39 
and Tf — 0.5. 

VI. APPLICATION TO NON-HOMOGENEOUS DNA 

As a last application of our method we shall consider the equilibrium configurations of non-homogeneous rods. We 
restrict ourselves to the simplest case of periodic non-homogeneities in Young's modulus. The motivation for this 
study is the fact that repeated (and therefore periodic) DNA sequences form a substantial fraction of all eukaryotic 
genomes |3^, The calculations presented here are based on the stiffness parameters recently computed for the 
32 tri-nucleotide units from DNA data | |4l[ |. Our goal is to understand how much the equilibrium configuration 
of a non-homogeneous rod differs from that of the homogeneous case when the rod is subject to fixed mechanical 
conditions. 

Repetitive DNA is formed by nucleotide sequences of varying lengths and compositions. Repeated sequences, 
reaching up to 100 megabasepairs of length | |40[ |, appear to have little or no functional role, and are commonly 
regarded as "selfish" or "junk" DNA |Q. We shall use a simple periodic formula for the (scaled) Young's modulus 
that covers most of the parameter interval spanned by the tri- nucleotides given in ref . p| : 

E{s) = l + acos—s . (52) 

£ is the period of the oscillations of the Young's modulus along the DNA and a is a parameter that depends on the 
specific sequence being repeated and that can not be greater than 0.66. 
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Eqs.(p2|)-(p5|) have to be slightly modified to include the non-constant Young's and shear moduli. We obtain 




(53) 



with 



Pe 



E{s)e', 



(54) 



Af3 = roM(s) ((/)' + 7/-' cos 0) 



(55) 



= Mz = E{s)^' sin^ e + P^ cos 6*. 



(56) 



The solutions do not depend on /i(s), since it does not enter in the differential equations for 6, Pg and ?/;■ Notice 
that these equations are not integrable if a 7^ 0. Although P^ and are still constants, the elastic energy per unit 
arclength is not. 

The method for solving the BVP for a non-homogeneous rod is the following: consider a solution extending from 
the origin to r/ with a = and initial conditions 6^ and Pg^ . Now integrate the Kirchhoff equations from the same 
initial conditions but using Eq.(|5^) with a 7^ 0. This new solution, whose end point is r/ -I- 5v, can be used as a trial 
solution for the non-homogeneous rod. Using the method of section IV we push the rod back to r/ . 

Figure || shows the up and down solutions for load parameters P = 1, P^ = 0, P^ = 10, r/ = 0.9 and z/ = for 
the homogeneous rod (black curve) and a non-homogeneous rod with a — 0.66 and C — 0.1 (gray curve). 

Finally, Figure ^ shows the effect of changing the period of oscillations in Young's modulus. We show the up and 
down solutions for non-homogeneous rods with the same load parameters of the previous figure and a — 0.66. The 
curves show rods for £ = 0.1 (gray), £ = 0.5 (thick black) and £ — 0.65 (thin black). The changes in the three- 
dimensional shape of the rods are evident for these values of load parameters. The sensitivity of the three-dimensional 
shape to the base-pair sequences indicated that DNA repeats may have conformational roles. 



In this work we presented a general method to solve the boundary value problem (BVP) for finite Kirchhoff 
filaments. The method consists in making small changes to a known trial rod that satisfies the Kirchhoff equations 
but not necessarily the boundary conditions. We combine a shooting technique with the method of monodromy matrix 
to push the end point of the trial solution to the desired position, step by step. By linearizing the Kirchhoff equations 
we obtain an explicit relation between a variation of the initial conditions (expressed in terms of Euler angles) and 
the consequent variation of the rod's end point. 

The solutions of the BVP are limited by the physical constraints of the rod, such as the moments and tension. 
A sketch of the allowed end points can be constructed by integrating the Kirchhoff equations for a large number of 
initial conditions. The regions of possible end points form complex figures reflecting the non- linear character of the 
equations. The regions of existence of end points may serve as a guide to find the appropriate load parameters needed 
for a desired solution. 

We presented several examples of rods with different end positions at different distances from the origin for various 
sets of load parameters. In all cases the method worked very well and the BVP was solved with at least six significant 
figures in the values of rf and Zf. We also applied the method to non-homogeneous, sequence-dependent, DNAs. 
We modeled pieces of repeated sequences by a sinusoidal oscillation of the Young's modulus. We showed that the 
tri-dimensional structure of the DNA is indeed sensitive to the presence of such sequences, a property that has been 
considered before |Q but studied only in terms of the DNA intrinsic curvature [Q. The effect studied here may 
contribute to other sequence-dependent properties that affect the three-dimensional conformations of the DNA. 



VII. CONCLUSIONS 
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FIG. 1: Trial (thin line), intermediate (thick line) and final (thick gray line) solutions for a rod with desired end point at 
r/ = (0.7, 0.0). The load parameters are = 1, P4, = 1 and F = 1, in scaled units. 



1.00 




FIG. 2: Regions of the existence of final points for = 0, = 1, F ^ 0.1. The curves D = 1.0 (fuU thick line), D = 0.9 (full 
line), D = 0.7 (dashed line) and D = 0.4 (dotted-dashed line) are also shown. The circles correspond to the rods in the Figure 
|[The plot was generated by varying the initial conditions in the intervals < 9o < n (30 points) and —5.0 < Peo < 5.0 (200 
points). 
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FIG. 3: Up and down solutions (drawn in the same color) for the end points shown by the circles in Fig. g. (a) D — 0.4; (b) 
D = 0.7; and (c) D = 0.9. 
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FIG. 4: Regions of the existence of final points for = 0, = 5, F — 1. The curves D = 1.0 (full thick line), D = 0.9 (full 
line), D — 0.7 (dashed line) and D — 0.4 (dotted-dashed line) are also shown. The circles indicate the rods drawn in Figure 
The plot was generated by varying the initial conditions in the intervals < 9o < n (30 points) and —20.0 < Peo < 20.0 (200 
points). 
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FIG. 5: Rods for = and F — 1. The lines corresponds to = 1 (thin black), P^ = 5 (thick black) and P^ — 10 (gray), 
(a) up solutions for r/ = 0.7; (b) down solutions for r/ = 0.7; (c) up solutions for r/ = 0.9; (d) down solutions for r/ = 0.9. 
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FIG. 6: Regions of the existence of final points for = = 5, F = 1. The curves D = 1.0 (full thick line), D = 0.9 (full 
line), D = 0.7 (dashed line) and D = 0.4 (dotted-dashed line) are also shown. The circles indicate the rods drawn in Figure ^. 
The plot was generated by varying the initial conditions in the intervals < 9o < n (30 points) and —10.0 < Peo < 10.0 (200 
points). 
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FIG. 7: Rods for P-i/j = and F = 1. The lines corresponds to P4, = P4, = 10 (black) and P4, = P4, = 5 (gray), (a) and (b) 
show the up and down solutions for r/ = 0.39 and Zf = 0.8; (c) and (d) show the up and down solutions for r/ = 0.5 and 
Zf = 0.49 
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FIG. 8: Comparison between homogeneous and non-homogeneous rods. Panels (a) and (b) show the up and down solutions 
respectively. The parameters are = 0, P^, = 10, F = 1, r/ = 0.9 and Zf — 0. The black curve shows the homogeneous rod 
and the gray curve the non-homogeneous rod with a — 0.66 and £ = 0.1. 




FIG. 9: Comparison between non- homogeneous rods with different periods C. Panels (a) and (b) show the up and down 
solutions respectively. The load parameters and final position r/ are the same as in Fig. @. The curves show rods for £ = 0.1 
(gray), £ = 0.5 (thick black) and £ — 0.65 (thin black) 



